Thermal diffusion of supersonic solitons in an anharmonic chain of atoms 
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We study the non-equilibrium diffusion dynamics of supersonic lattice solitons in a classical chain 
of atoms with nearest-neighbor interactions coupled to a heat bath. As a specific example we choose 
an interaction with cubic anharmonicity. The coupling between the system and a thermal bath with 
a given temperature is made by adding noise, delta-correlated in time and space, and damping 
to the set of discrete equations of motion. Working in the continuum limit and changing to the 
sound velocity frame we derive a Korteweg-de Vries equation with noise and damping. We apply 
a collective coordinate approach which yields two stochastic ODEs which are solved approximately 
by a perturbation analysis. This finally yields analytical expressions for the variances of the soliton 
position and velocity. We perform Langevin dynamics simulations for the original discrete system 
which confirm the predictions of our analytical calculations, namely noise-induced superdiffusive 
behavior which scales with the temperature and depends strongly on the initial soliton velocity. 
A normal diffusion behavior is observed for solitons with very low energy where the noise-induced 
phonons also make a significant contribution to the soliton diffusion. 

PACS numbers: 05.40.-a,63.10.+a,05.45.Yv 



I. INTRODUCTION 



Nonlinear onc-dimcnsional lattice dynamics, namely 
propagation of coherent excitations in monatomic chains 
modeling discrete microscopic structures, is associated 
with several important problems in physics. Among these 
excitations are solitary waves, which for simplicity are 
called here solitons. These solitons can be supported 
by chains with realistic interaction potentials between 
the particles |. They are supersonic non-topological 
collective excitations. In spite of their relative simplic- 
ity, the solitons clarify many features of molecular chains 
p, |J. ||, [| 0] . For example, due to their robust character, 
lattice solitons have been used to model the energy trans- 
port in polypeptide chains in muscle proteins |], |^, [n]] 
or the energy transport in DNA p"2[ . Numerical simula- 
tions at realistic temperatures for transport in proteins 
have shown that lattice solitons can propagate over long 
distances in a chain with the Lennard- Jones potential Q . 
Moreover, the lattice solitons are more stable than Davy- 
dov solitons if collisions between the two types of solitons 
are considered . There is no clear evidence that lattice 
solitons like a Toda type, which are non-topological, can 
exist in thermal equilibrium. This holds both for static 
properties, like the specific heat, and for dynamics quan- 
tities, like the dynamic form factor (Fourier Transform 
of the displacement autocorrelation) . On the other 
hand, there exists evidence from real experiments that 
strain solitons can be generated and observed in non- 
linear elastic rods flial. These solitons in some cases can 



be described by Korteweg-de Vries (KdV) type solitons, 
which are non-topological. 

To our knowledge there are no previous analyti- 
cal studies supported by Langevin simulations about 
non-topological lattice soliton diffusion in anharmonic 
monatomic chains of particles with nearest-neighbor in- 
teractions. There are many studies on stochastic partial 
differential equations, in particular stochastic KdV-type 
equations have been extensively studied numerically and 
analytically (l4[ [TJ, [l6], [l7| |8| |l9) due to the integrabil- 
ity of the KdV equation. In fact, the KdV equation is a 
good approximation to describe analytically the dynam- 
ics of lattice solitons on a monatomic chain with nearest- 
neighbor interaction and cubic anharmonicity if the soli- 
ton velocity is very close to the sound velocity (very-low- 
energy solitons) j|, . Notice that for a polynomial po- 
tential, namely harmonic term plus cubic or/and quartic 
anharmonicity, the one-soliton solution of the KdV equa- 
tion is known analytically, while for more realistic inter- 
action potentials like Lennard- Jones or Morse there are 
no exact soliton solutions. In the more general context of 
lattice systems, there are a few analytical studies about 
diffusion of coherent lattice excitations, viz. stochastic 
vortex dynamics in two-dimensional easy-plane ferromag- 
nets J2l| or soliton diffusion on the classical, isotropic 
Heisenberg chain |22l . 

The aim of this work is to provide an approximate 
analytical description of the soliton diffusion dynamics 
in a monatomic chain with a cubic anharmonicity under 
thermal fluctuations. For this purpose we generate a sin- 
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gle soliton which has an energy far greater than kg T, 
where kg is the Boltzmann constant, and T is the tem- 
perature of a thermal bath. This soliton propagating on 
a chain in contact with the thermal bath shows a dif- 
fusive behavior. We consider this soliton diffusion dy- 
namics during the thcrmalization process of the system 
up to times when the system energy has relaxed nearly 
to its stationary value. This means that we study the 
non-equilibrium diffusion dynamics of lattice solitons on 
anharmonic chains subject to thermal fluctuations. 

In order to perform the coupling between the system 
and a thermal bath with a given temperature an addi- 
tive noise term, providing energy input, is added to the 
discrete equations of motion. This term has to be bal- 
anced by a damping term providing energy dissipation. 
Here, we suggest as a damping term the so-called hy- 
drodynamical damping pjj which is extensively used in, 
e.g., elasticity theory. Notice that this type of damping 
is due to irreversible processes taking place within the 
system. The corresponding noise term, which fulfills the 
fluctuation-dissipation theorem, takes the form of a dis- 
crete gradient of Gaussian white noise delta-correlated in 
space and time. A similar Langevin-type equations has 
previously been considered in the context of mesoscopic 
Langevin dynamics p5[ . 

We notice that our system in the continuum limit can 
be approximated by a noisy KdV-Burgers-type equation 
pG| , B7| . So in this case we can use the one-soliton solu- 
tion of the KdV equation not only as initial condition of 
our discrete system but also in our analytical approach 
in the continuum limit. Notice that the shape of broad 
KdV solitons tends to be identical to the shape of broad 
supersonic lattice solitons 0. In this work we apply a 
generalized traveling wave ansatz combined with a collec- 
tive coordinate formalism in the framework of the KdV 
equation as an analytical approach to study the diffusion 
of lattice solitons. 

In the next section we present the equations of mo- 
tion of our discrete system. From this we formulate a set 
of stochastic equations of motion by adding noise and 
damping. Next, we apply the continuum limit and de- 



rive a form of noisy KdV-Burgers equation. In section 
we apply a collective coordinate approach which yie 



|ffl 

ds 



analytical expressions for the thermal averages and vari 



ances of the soliton position and velocity. In section IV 



we compare our analytical predictions with the results 
from Langevin dynamics simulations for the original dis- 
crete system. Our conclusions are summarized in the last 
section. 



II. THE CONTINUUM LIMIT 

We consider an anharmonic chain of particles with 
mass M and nearest-neighbor interactions. The parti- 
cles interact via an anharmonic potential with a cubic 



anharmonicity. The Hamiltonian of this system reads 
[Pi 



H 



E 



\2M 



G ( \{Yn+i - Y n f + j(Y n+1 - Y n f 



(1) 



where Y n denotes the longitudinal displacement of the 
n-th particle from its equilibrium position, and 



Pn=M 



dYn 

dt 



(2) 



is the momentum. Here G and A are the potential param- 
eters whose values depend on the lattice. The associated 
first order equations of motion read 



dY n 
dt 

dP n 

dt 



1 



:Pn 



dH 

dYn 



-^Damping 



(3) 
(4) 



where 



— = -G(Y n+1 -2Y n + Y n _ l ) 
oY n 

~GA ((Y n+1 - Y n ) 2 - (Y n - y„_ x ) 2 ) 



(5) 



In Eq. we have already added both a stochastic force, 
pNoise^ anc j a damping force, F® ampmg , Both forces cou- 
ple the discrete system with a thermal bath. Here, we 
use the inner or hydrodynamical damping, which reads 

MM 



pDamping _ 



fdY r 



n+1 



dY n dY n -i 



\ dt 



dt 



dt 



(6) 



This means that the energy dissipation is provided by 
the irreversible processes arising from the finite velocity 
of the internal motions of the system, namely time deriva- 
tive of the relative displacements between particles in the 
chain. Eq. (^|) is the discrete version of the damping 
used in elasticity theory p4fl . To fulfill the fluctuation- 
dissipation theorem the noise must have the form (see 
App. A) 



F 



Noise 



where 



D (£„+i(t)- £„(*)) 



D = 2Mvk B T 



(7) 



(8) 



is the diffusion constant and v is the damping constant. 
£ n (t) is delta-correlated white noise, 



S nm S(t -t), 
0. 



(9) 
(10) 



Since our interest is the study of the lattice soliton dif- 
fusion close to the sound velocity, c, we can use the 
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continuum limit approach, where Y n (t) — > y(x,t) and 
£n(t) ~ * with x = na and a the equilibrium atomic 

spacing. In this limit Q, Eq. (0) reduces to a form of 
noisy and damped KdV equation (see App. B for details) 



d T u + 6ud s u + dgii = h>id ss u — \/Did s £(s,T) (11) 
where 

s = a(x — ct), t = f3t, u = jd s y. (12) 



The constants a, (3 an d 7 ar e defined in (B9), and v\ and 
Di are given by Eqs. ( B13[ ). Note that 



(as,TW,T')) = 6(s-s')5(T-r'). 



(13) 



Here and in the following the line over £ is omitted. 

The case D\ = reduces Eq. ( pi] ) to the KdV-Burgers 
equation. The associated KdV equation is 

d T u + 6ud s u + d 3 s u = (14) 

whose one-soliton solution reads 

u (s,t) = 2^ S ech 2 [77o(s - irfir - s )}. (15) 

Here 

Vo = -\/3c(v - c) (16) 
P 

is the inverse soliton width and so is the initial soliton 
position. The sound velocity c and the constant p are 
defined in (B5). 



III. COLLECTIVE COORDINATE APPROACH 

To analyze our problem we assume that the soliton pro- 
file, uq(s, t), is not disturbed by the noise and damping 
terms and that only the width and amplitude are mod- 
ified. This assumption is well satisfied for low-energy 
solitons, whose velocity is close to the sound velocity, be- 
cause tails induced by the perturbations are small in this 
velocity regime |2?| |29| . So we introduce a generalized 
traveling wave ansatz of the form 



u(s,t) = u Q (s - S(r),ri(r)) 

= 2ry 2 (r) sech 2 [r/(T)(s - S(t))}, 



(17) 



where the collective variables S(t) and t](t) are the soli- 
ton position and the inverse soliton width, respectively. 
Here and in the following the index of the one-soliton 
solution uq is omitted. 

To obtain the equations for our collective coordinates 
we follow Q |3|. First, by substituting @ into Eq. 
(0) we get 

4>i S{t) + 4> 2 r)(r) = v x d ss u - y^xd^s, r), (18) 



where 



and 



<Ms,t) 



<Ms,t) 



t)ii 
dS 



du 
dr]' 



(19) 



(20) 



Notice that the functions {4>i\ i= i 2 coincide with the adi- 
abatic approximation (omitting secular terms in time) 
of the discrete solutions of the linearized KdV equation 
around the one-soliton solution (|l^) |32]]. We remark 
here that our collective coordinate theory does not take 
into account the contribution of the phonon modes (con- 
tinuous basis function solution of the linearized KdV). 
We d iscuss the effect of noise-induced phonons in section 
jVB. The functions {4>i} i=1 2 are a l so orthogonal, so the 
inner product J ds (/>i(s,r) h(s,r) projects a function h 
onto the functions {4>i} i=1 2 - Then, by projecting Eq. 
@ we get 



A i S(T)+B i f,(r) = f i + ft ampin9 - 
where 



7, 



Nc 



= / ds 



A, 
B, 
fi 



rdampinq 



p Noise 



_ du 



du 



i = l,2 , 
(21) 



(22) 
(23) 



07] 

ds (6ud s u + d^u) fa, (24) 



dsd ss u 4>i 



D 1 / dsd s £(s,T) 



(25) 
(26) 



After some calculations the Eqs. (21) take the form 



dS{r) 
dr 

dr 



4t7 2 (t) + 



15VA 



30^i 



30 + 7T 

45A/DT 



64?7 5 (r) 



ds(ds<l>iMs,T), (27) 



16(30 + tt 2 )77(t) 



(28) 



To achieve the calculations we have assumed that the soli- 
ton profile remains mostly unaffected and only its width 
and amplitude change due to the stochastic perturba- 
tions. Then, at least for small noise, we can perform the 
calculations by taking the soliton field out of the aver- 
ages. Moreover, we have interpreted Eqs. (|27]) and (|28| ) 
in the Stratonovich sense, because it assumes £(s,t) is 
a real noise with finite correlation time, which is then 
allowed to become infinitesimally small after calculat- 
ing measurable quantities [ p3| . Notice that white noise 
means taking the limit of zero correlation time. 
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Eqs. (E7J) and fcq) can take the form 



dY(r) 
dr 



= A ifr [Y(r)]+ d S B^[ S ,Y(r)]£( S ,r), (29) 



where the elements of the noise vector £ satisfy (13). 
{Yi, >2} = {S 7 rj} are the elements of the vector Y, the 
elements {Af tr , Af tr } of the drift vector A str are the 
drift terms in Eqs. (|27j) and (|2^), respectively. The dif- 
fusion matrix B s ' r is diagonal, where Bf* r and B^ r are 
the coefficients in front of the noise in Eqs. (27) and 
(p8|), respectively. In order to facilitate the calculations 
we write Eq. in the Ito form, 

dY(r) = A /to [Y(r)]dr+ / <2sB Ito [s, Y(r)]dW(s, r), 

(30) 

where the dW(s, r) = £(s, t)c£t is a Wiener process. Via 
a Fokker-Planck equation, one can show that the ele- 
ments of the drift vector A Ito read (3j| 

A[*°[Y(r)]=AfqY(r)] + 

~£ / d S S^[ S ,Y(r)]^Bff[ S ,Y(r)] 



while 



j km' 

hj,k= 1,2 



B J *°[ S ,Y(r)] =B^[ S ,Y(r)] 



(31) 



(32) 



Notice that A Ito and B J *° are nonanticipating functions. 
So, from Eq. (|30|) it is easy to show the following averages 



(S(r)) = ( J dr'W)} 



T 

, 30z^i 3 



fa(r)> 



dr' 



30 + 7T 2 



r/r 



,2 25(231 + 8tt 2 )£>i 
112(30 + 7T 2 ) 2 



Var(7?(r)) = ( / dr 



, 225(21 + 7t 2 )Di?7(t') 



28(30 + 7T 2 ) 



2\2 



Now we define a new set of Langevin equations, 

dY»(r) = a 4 dr + 6 y - dWj(r) 

with i,j = l,2 and {Y x , Y 2 } = {S, rj}, (35) 

which we have interpreted in the Ito sense. dWj (r) = 
(r)dr are Wiener processes where we have let the noises 
to be uncorrelated, namely 



Mr)^{r'))=5{r-r')8 nl . 



(36) 



In order to determine the values of a* and 6i, we have 



demanded that Eqs. (|35|) satisfy the relations (p3[). It is 
straightforward to see that Eqs. (|35| ) take the form 



dS(r) = 4r; 2 (T)dr + 
dr)(r) = 



5^3 / L»i 



4\/7 V »7 3 (t) 



dWi(r) 



30i/i 



30 + 7T 



225(231 + STr 2 )^! 
112(30 + 7T 2 ) 2 



15V2T 



2\/7(30 + Tr 2 ) 



(37) 



dr 



(38) 



Eqs. ( |3Th a nd fl38|) are statistically equivalent to Eqs 
( p7| ) and (|28|) because they share the same Fokker-Planck 
equation. Since the derivation of Eqs ( |37j ) and (^8|) in- 
volved approximations, we have not solved them exactly. 
Instead of that, we have used perturbation analysis [p4[ , 
which is developed in detail in App. C. In order to do so, 
we have considered the thermal terms as perturbations, 
so Eqs. d37|) and (Esb take the form 



dS(r) = A V 2 (r)dr + e^L/ ^TT dWi(r) 



4\/7 V ?? 3 (t") 



dq{r) 



30^i . , 

3oT^ (r)dr + e 



/225(231 + 8tt 2 ) j Di 



V 112(30 + 7T 2 ) 2 



15V21 + 7T 2 
2^7(30 + Tr 2 ) 



(39) 



dr 



(40) 



Now, we seek an asymptotic solution in the form of a 
small-noise expansion 



S(t) = s (r) + e si(t) + ■ 
= %(t) +er?i(r) + • 



(41) 



CW(S*(t)?7(t)) = 0. (33) 
Here ( ••• ) means average over an ensemble of realizations, 



Corr(PQ) = (PQ)-(P)(Q) and Var{P) = Corr(PP) 

(34) 



Here, e is a factor introduced for convenience in the an- 
alytical calculations. Notice that the case e = reduces 
Eqs. (|39|) and (Efl) to the deterministic case. In order 
to interpret our perturbation theory we must set e = 1 
and assume that the terms on the r.h.s. of Eqs. ( |39| ) and 
( (40| ) are sufficiently small. So we must restrict ourselves 
to a regime of low temperatures of the thermal bath (Di 
small). From the perturbation analysis we obtain that 
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(S(r)) = ( So (r)) + ( S i(r)) 

, % 2 (0) , n , , , , 15-Pi (231 + 8^ 2 ) 7ft,(0) (2(1 + Ar) 5 / 2 - 5Ar - 2) 
= 4— log (1 + At) + , 

(A V 2 (t)) = (4 rf (r) + 8tjo(t)»ji (r)) 

4?7 2 (0 ) 45^ 1 (231 + 87r 2 )r ? o(0)((l + Ar) 5 / 2 - 1) 



Var(S(T)) = £>i 



Var(4 r? 2 (r)) 



1 + Ar 7(30 + tt 2 ) 2 A(1 + Ar) 2 

-15 480 (21 + 7r 2 )r7g(0)(8 + 7At(5At + 4)) 



56 ,7 3 (0) A 49 (30 + tt 2 ) 2 A 3 (1 + A t) 2 

15(1 + T 2 A)3 , /2 (2048(21 + 7r 2 )ry 6 (0) + 7(30 + tt 2 ) 2 A 2 (1 + Ar)) ) , 

392 (30 + 7r 2 ) t)q(0) A 3 V ° V 



7200L>i(21 + 7r 2 )?7 3 (0) 



49A(30 + tt 2 ) 2 Vv/TTA^ (1 + Ar) 4 



(42) 



The expr essions for D\ and A are given by Eqs. ( |B13| ) 
and ( C1C| ), respectively. 



IV. SIMULATIONS 

Substituting Eqs. (^]) and (0) into Eq. (||) we get the 
full set of discrete equations of motion written in abso- 
lute displacements. However, for our simulations relative 
displacements are more convenient, because the lattice 
solitons in this representation are pulse solitons whose 
amplitude vanishes at infinity. This characteristic allows 
us to use periodic boundary conditions which are neces- 
sary for long simulation times, because we want to avoid 
reflections at the boundaries. So the discrete equations 
of motion in relative displacements read 



M 



dt 2 



G(V n+1 -2V n + V n - 1 ) + 



GA (V, 



71+1 



2V 2 + V 1 

M ' n ' v n—l 

dV n dV, 



> r ( dVn+i n av n uv n -i 

Mu — 2— — 

V dt dt dt 



D(£ n+ i(t)-2Ut)+tn-i(t)), (43) 

where V n (t) = Y n+1 (t) - Y n (t) and D = 2Mvk B T. The 
periodic boundary conditions read 



d l V S V N -x 



d l V N d l Vi 



dt 1 



dt 1 



dt 1 



dt 1 



1 = 0,1 



£oW = 6v-iW, 6v(t) = &(i), 



(44) 



where N is the number of particles of our chain and N—l 
is the number of bonds. 

A suitable method to detect the position of a pulse 
lattice soliton, V n (t), is to search for its maximum p8|. 



However, in the presence of stochastic perturbations this 
method is not useful since the pulse shape is strongly 
masked by the noise, an example of this situation is 
shown in App. So from the data of our simulations 
we have taken snapshots of the system at different times, 
and from them we have generated the kink shape Y n {t) 
of the lattice soliton by using the algorithm 



Y n {t)=Y 1 {t) + Y,V i {t) n = 2,3,---,N. (45) 

i=l 

The kink shape is less distorted by the noise than the 
pulse shape Vi(t). In ( (45| ) Yi(t) is a boundary condition 
that we have demanded to be 



N-l 



(46) 



so at t = the amplitude of the center of the kink shape 
is zero, as it should be from the theory B. Notice that 



N-l 



Y N (t)-Y 1 (t) = J2 V ^ 



(47) 



is a conserved quantity in our system, i.e. 

N-l 



Y N (t)-Y 1 {t)= J2^(t) = 0. 



(48) 



i=l 



We have checked Eq. ([48]) with a precision higher than 
10 -14 over the whole time range of our Langevin dynam- 
ics simulations. 

In order to determine every time the parameters of the 
soliton, namely soliton velocity v and position x, we have 
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proceeded as follows. We have searched for the values of 
the parameters x and v where the relation 

no+rii 

(Y n -y (na-x,v))=0 (49) 

n—no—ni 

is fulfilled. Here 



. . 6\/2 h civ — c) / na — x , 
yo{na — x,v) = tanh 1 (50) 



L{v) 



with 



L(v) = (a^t))- 1 



2 c (v - c) 



Here the function n is def ined in ( lq) a nd the constants h, 
p, c and a are defined in (B5) and ( p9| ). the function j50| ) 
is the absolute displacement representation of the one- 
soliton solution ( |15|) in a frame moving with the soliton 
velocity. In Eq. (E9j) 



no = int(x) and n\ 



int(-L( V (0))), 



(52) 



where int(-) denotes the integer part of a number and v(0) 
is the initial soliton velocity. The value of n\ has been 
chosen to take into account only the core of the lattice 
kink-shape and it is constant during our simulations. In 
order to determine both x and v we consider different 
values of v in Eq. ( ^9] ) within a range of velocities around 
the initial soliton velocity, namely v — c 6 [0.1 (i>(0) — 
c), 2 (v(0) — c)]. For every value of v we search the value 
x that fulfills Eq. (j49|), so we get a set of pairs of values a; 
and u. Finally, from this set of pairs of values we search, 
by using linear interpolation, the values of x and v which 
fulfill the relation 



no+ni 

Y n y (na- x,v) 

n=riQ—ni 
na+ni 

J2 (yo(na- x,v)) 2 

n—riQ — n i 



(53) 



Notice that in Eqs. (^9|) and (|53|) we have assumed that 
the lattice kink shape, Y n , is closely related with the 
function y , however, as was mentioned in Ref. p^ j, a 
pulse lattice soliton in the presence of damping develops 
a tail. The amplitude of this trailing tail depends on both 
soliton the velocity and the damping, so it is bigger when 
the damping and/or the soliton velocity is higher. Thus, 
we restrict ourselves to velocities very close to the sound 
velocity where the effect of this trailing tail is negligible. 

Up to now we have determined the parameters x and 
v, which fit the function j/o to the lattice kink shape Y n , 
so we have not measured directly either x or v. Since the 
function yo is closely related to the lattice kink-shape 
Y n , one could assume both x and v as an estimate of the 
soliton position and velocity, respectively. However, we 
have taken only the parameter v as an estimate of the 



soliton velocity and with this value we have used a dif- 
ferent method to determine the soliton position. In fact, 
in order to be in agreement with our collective coordi- 
nate approach, where we have projected the equations 
of motion onto the Goldstone mode <j>\ (Eq. ([l9|)), we 
have projected the noisy kink shape Y n (t) onto the pulse 
solution wo defined in (|l5|). Notice that in the absolute 
displacement representation the function uq is the Gold- 
stone mode. So this projection reads 



n+n 2 

P(x) = ^ Yi(t)u (ia- x,v) 



(51) 

where 



uo(i a - 



6c(v 



— sech 2 



P 



i a — x 



(54) 



(55) 



and x — na. The value of n-i in Eq. ( |54| ) is much larger 
than the soliton width, so the boundary effects are neg- 
ligible. The function uq{i a — x, v) is the one-soliton so- 
lution ( |l5| ) in a frame moving with the soliton velocity. 
Afterwards we have searched, by linear interpolation, the 
value x where P(x) vanishes and we have defined it as 
the position of the soliton center of mass. At this point 
we remark that the values of x following from this latter 
method are not significantly different from the values of 
x following from the former method (Eqs. ^ and [33]). 
However, we consider the latter method to be more ap- 
propriate than the former one in the sense that we pro- 
ceed in our code in a similar way as in our analytical 
calculations. 

Our Langevin dynamics simulations were performed 
for a chain with 1500 lattice points. The time integration 
was carried out by using the Heun method |35| , which has 
been successfully used in the numerical solution of par- 
tial differential equations and difference-differential equa- 
tions, coupled to either an additive or a multiplicative 
noise term [^l], [22| p3[ pq| . Here, we have used the con- 
served quantity (p:7|jto check the accuracy of our code 
p8[ . For the longest simulation time the variation of this 
conserved quantity has been lower than 10 -9 %. In or- 
der to start the simulations at t = we have used the 
one-soliton solution (JlJ) of the KdV equation in the lab- 
oratory frame. The average values have been calculated 
over 200 realizations up to a final time 5000. All values of 
the constants of the equation ( [43| ) are set at unity except 
the damping constant which is set at v = 0.003. Notice 
that for lower values of damping the relaxation of the 
system energy would take more time in our simulations 
to reach a regime close to its stationary value. On the 
other hand, higher values of damping can strongly distort 
the soliton shape, namely the tail induced by the damp- 
ing cannot be neglected when the value of the damping 
is high. In App. |D| we show the thermalization process 
in our system. In our simulations the values of tempera- 
ture, T, and initial soliton velocity, w(0), are parameters 
(see figure captions). 
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FIG. 1: Averaged soliton position (a) and velocity (b) vs. 
time in the sound velocity frame z = x — t (c = 1), with 
v = 0.003 and T = 5 x 10 -5 . Dotted lines: simulation, solid 
lines: theory (Eqs.(^)). A, B, C, D, and E correspond to 
different initial velocities, namely v(0) = 1.001, 1.003, 1.005 
and 1.007, respectively. 



A. Soliton propagation 

In Figs. [l]a and ^> we show several examples of both 
the averaged soliton position, (x(t)}, and the averaged 
soliton velocity, (v(t)}, as functions of time from both the 
simulation and the theory (see Eqs.(M2)). Notice that 



and 



( x (t)) = -(S([3t))+ct 
a 



(56) 



(57) 



For all cases the soliton position from the simulation 
agrees well with the position given by the analytical the- 
ory. In the case of the soliton velocity, the agreement 
is better for initial velocities close to the sound veloc- 
ity than for higher velocities. In fact, the time evolution 
of the velocity given by the simulation is always higher 
than the theoretical prediction. This systematic differ- 
ence may be due to the small tail that is generated by 
the lattice soliton since it is a non-topological soliton p8| . 
So it may affect our numerical method for determining 
the soliton position. We point out that the amplitude 



of this tail depends on both the soliton velocity and the 
damping constant, and can be neglected only for veloc- 
ities close to the sound velocity and small values of the 
damping constant [ ^8| . This is the most important rea- 
son for restricting our study to low-energy solitons whose 
velocities are close to the sound velocity. Since this dif- 
ference is systematic it does not play any role in the nu- 
merical calculation of the variances which is our more 
important goal in this article. 



B. Soliton diffusion 

In Figs. H and || we show the variances of the soliton 
position and velocity vs. time for different initial veloci- 
ties. The temperature of the thermal bath in Fig. || is 10 
times higher than that in Fig. ^. The results scale very 
well by a factor of 10; i.e. the variances are pro portional 
to the temperature, as expected from Eqs. (|42|). 

Notice that our theory (solid lines) has no adjustable 
parameters. Taking into account that the theory consists 
of several steps (discrete system — > Bq equation — + KdV 
equation — > collective coordinate theory — > perturbation 
analysis) it is already a significant success to obtain the 
observed results, as seen in Figs. || and || 

We observe that the behavior of the variances depends 
strongly on the initial soliton velocity. For low-energy 
solitons, whose velocities are close to the sound velocity 
(Figs. I and | cases (a) and (b)), the soliton diffusion 
tends to be nearly normal, i.e. linear in time. In fact, 
our theory predicts a normal diffusion for times 



t « t* 



^30 



495 



v(v(0) - c) 



(58) 



This estimate was obtained by comparing the first with 
second terms of the Taylor expansi on in powers of r of the 
variance of the velocity (see Eq. ( C16[ )). For low-energy 
solitons (u(0)^ c) t* is much larger than our simulation 
time (Figs, ya and ^a). This means that the superdiffu- 
sivity is not very dominant. However, for higher-energy 
solitons the anomalous behavior turns out to be impor- 
tant after some time. In those cases t* is comparable 
with our simulation time (Figs. |2|e and ||e). 

In the case of solitons with very low energy the vari- 
ance of the position (Figs. |^a and ||a) does not agree so 
well with the theoretical prediction. In fact, we observe 
a transient behavior for times t<Z/v = 1000 where the 
system energy shows a fast relaxation process (see App. 
Pi) . Those discrepancies between theory and simulations 
may be due to the combination of two effects. First, the 
profile of low-energy solitons is strongly masked by the 
noise, so the numerical detection of the position can be 
distorted. Second, not only the noise but also the noise- 
induced phonons can make a significant contribution to 
the variance of the position, since the reduced tempera- 
ture, T = Ub T/H(0) (temperature in units of the initial 
soliton energy H(0)), of the thermal bath is higher here 
than in the other cases (see captions Figs. || and ||) . In 
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FIG. 2: Variances of the soliton position (panels a, c and e) and velocity (pane ls b, d an d f) o f the soliton vs. time, with 
v = 0.003 and T = 5 x 10~ 6 . Dotted lines: simulation, solid lines: theory (Eqs.(C13) and (|C14j)). The panels correspond to 
different initial velocities, namely v(0) = 1.003 (a and b), 1.005 (c and d) and 1.007 (e and f). The reduced temperatures are 
T = 0.0061419, 0.00283324, 0.00169765, respectively. 



this respect we estimate the phonon effect on the diffu- 
sion of low-energy solitons in the next section (see also 
Fig. @). 

Notice that the variance of the soliton position is larger 
for low-energy solitons (Figs. ||a and||a) than for higher- 
energy solitons (for instance Figs. ||s and^e). This is due 
to the fact that the higher-energy solitons are more ro- 
bust against thermal fluctuations than the lower-energy 
ones. Or, equivalently, the reduced temperature T of the 
thermal bath is higher for slow solitons than for the fast 
ones (see captions of Figs. || and ||). 

On the other hand, the superdiffusive behavior is more 
pronounced for higher-energy solitons. This is because 
the soliton velocity turns out to be more sensitive to 



the thermal fluctuations in this case than in the case of 
broader solitons. Notice that the soliton velocity and soli- 
ton width are related. Also, since the higher-energy soli- 
tons encompass few lattice sites, the soliton- width per- 
turbations are larger with respect to the averaged soliton 
width in this case than in the case of broader solitons 
(low-energy solitons). In fact, the variance of the soli- 
ton velocity shows this effect, namely that for broader 
solitons (Figs. J^Jd and ||b) this variance is smaller than 
for narrower solitons (for instance Figs. ||f and ^f). The 
discrepancy between our theory and the numerical sim- 
ulations for higher-energy solitons (Figs. and ||, cases 
(c)-(f)) is mainly due to the fact that our theory is valid 
only for soliton velocities close to the sound velocity. 
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FIG. 3: Variances of the soliton position (panels a, c and e) and velocity (pane ls b, d an d f) o f the soliton vs. time, with 
v = 0.003 and T = 5 x 10 -5 . Dotted lines: simulation, solid lines: theory (Eqs.(C13) and (|C14j)). The panels correspond to 
different initial velocities, namely v(0) = 1.003 (a and b), 1.005 (c and d) and 1.007 (e and f). The reduced temperatures are 
T = 0.061419, 0.0283324, 0.0169765, respectively. 



With respect to the variance of the soliton velocity 
(Figs. H and |L panels b, d and e) we observe that it 
is mostly anomalous and its behavior is nearly quanti- 
tatively predicted by our theory for < t < 2000 in all 
the cases. For larger times, 2000 < t < 5000, there is a 
discrepancy which becomes larger with increasing initial 
velocities. 

We remark that the numerical results shown in Figs. ^ 
and H do not change for systems with the double number 
of sites, namely 3000. 

We comment that there was a previous attempt by 
Scalerandi et al. [[l9] to calculate theoretically the mean 
square displacement of a KdV soliton subject to stochas- 
tic fluctuations. They considered the case of small Stokes 



damping and a simple white noise delta-correlated in 
time and space. Though their theoretical result shows 
the appearance of a noise-induced superdiffusive behav- 
ior, it does not have the same dependence with respect 
to the soliton width as our results (^) , which agree well 
with our simulations. 



C. Estimate of the phonon contribution 

In order to estimate this contribution we have per- 
formed the following numerical test. We have simulated 
the propagation of low-energy solitons under thermal 
fluctuations up to a time (e.g. 2500) when the system 
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energy is close to its stationary value (see App. ^). No- 
tice that the diffusion of low-energy solitons is mostly 
normal (see Figs ||a and ||a). Afterwards, we have iso- 
lated the system from the thermal bath by switching off 
noise and damping, so that solitons propagate only in 
the noise-induced phonon bath. The diffusion is mostly 
normal before and after the system is isolated, i.e. the 
variance of the position is linear in time. We have com- 
pared the slope of the variance of the isolated system 
(t > 2500) with the slope of the variance in the case 
when the system is in contact with the thermal bath the 
whole time. An example of this test is shown in Fig. ||. 
We observe that both slopes, after switching off noise and 
damping, are different. In a normal diffusion process the 
slope of the variance of the soliton position is the diffu- 
sion constant. We call this Dtotai when the system is in 
contact with the thermal bath since there is a contribu- 
tion of both the thermal fluctuations and the phonons. 
The diffusion constant due to the noise-induced phonon 
bath (isolated system) is termed D p h- On the other hand, 
our theoretical diffusion constant, D^, is defined as the 
linear coef ficien t of t he T aylor expansion of Var{x(t)) 
(see Eqs. (C13) and (C14)). Since our theory does not 
take into account the contribution of phonon modes, we 
expect that the value D no i se = Dtotai — D p h may be 
of the same order of magnitude of D t h- We observe in 
Fig. ^| that the relative deviation of D t h from D no i se , 
{(D noise - D t h)/D noise , has the same order of magnitude 
of D t h which is not surprising since the soliton shape is 
strongly masked and distorted by the noise (see App. |]). 
In this respect we remark that from our results in Figs. 
H and H we observe that 



y/Var(x(t)) < L(v(t)), 



(59) 



where the soliton width L(v(t)) is defined in (J5l|). The 
relation (|5^) means that the stochastic deviations of the 
soliton center from its mean value are relatively small 
compared with the soliton width. So the diffusive dy- 
namics of the soliton position evolves inside the soliton 
core. Thus, this dynamics is very sensitive to the fluc- 
tuations of the soliton shape. Notice that our method 
of determining the soliton position depends implicitly on 
the soliton shape. So in the case of low-energy solitons, 
where the shape is strongly masked by the noise, the 
uncertainties of the method of soliton detection are rela- 
tively large compared with the diffusive dynamics of the 
soliton position. This is because the diffusive dynam- 
ics is relative small with respect to the soliton width. In 
this respect we note that higher-energy solitons present a 
well defined shape, i.e. the thermal fluctuations are small 
respect to the soliton amplitude, so the variance of the 
soliton position, given by our detection method, indeed 
agrees better with our theory, (see ||c and ||c). Finally, 
we stress that in our tests the phonon contribution to the 
soliton diffusion could be clearly observed only for very 
low-energy solitons (Figs. || and ^J, cases (c) and (e)), 
for higher-energy solitons the effect is negligible, namely 
D ph ~ 0. 




1000 2000 3000 4000 5000 



FIG. 4: An example of a test to determine the phonon con- 
tribution to the diffusion constant. Comparison of the be- 
havior of Var(x) from simulations (dotted lines) when the 
system is in contact (A) and isolated (B: noise and damping 
are switched off for t > 2500) from the thermal bath. The 
slope of the straight lines (solid lines) fitted to the simula- 
tion data (dotted) in both cases give the observable values of 
the diffusion constant, namely Dtotai (A) and D p h (B). The 
difference can be compared with the sl ope D t h of the dashed 
line (linear part of the theory, Eq. (|oil|)). v(0) = 1.003, 
T = 5 x 10" 6 , T = 0.0061419. 
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FIG. 5: Relative deviation ((D noise — D t h)/D noise ) of the 
diffusion coefficient vs. k B T/H(0) = T with v(0) = 1.003. 



D. Estimate of the bath temperature in real 
physical systems 

To clarify the physical meaning of the obtained results 
we have estimated the soliton energy and characteristic 
temperature of the thermal bath for two systems where 
solitons are believed to play an essential role: a-helical 
proteins || and crystal inertial gases [B7p. In both cases 
the Lennard- Jones (LJ) interaction potential, 



!>lj = 4£ ( ( - 



0' 



(60) 



is used. One can estimate the potential parameters G and 
A of the Hamiltonian ([!]) by using Taylor expansion of 
( |60| ) around the minimum r = (1 — 2 1 / 6 )a. So, taking into 
account only the coefficients of the harmonic and cubic 
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terms of this expansion, one can determine the potential 
parameters, namely 



A = - 



21 



and G = 



36 2 2 / 3 J B fl 



(61) 



By using the one-soliton solution of the Bq equation (Bl) 
Q and performing an integration instead of a summation 
in Eq. ([!]), the initial soliton energy H(0) reads 



H(0) = 



16V3 
245 



((«/c) 2 -l) 8/2 (l + 9(«/c) a )£b. (62) 



Here v and c are the soliton and sound velocities, respec- 
tively. £0 = 0.22 eV for a-helix § or E Q = 1-0 x 10" 2 eV 
for argon p8| . So, for example, in the cases of w/c=1.003 
or 1.03 we get 



v/c 


a-helix 
H(0)/k B [K] 


argon 
H(Q)/k B [K] 


1.003 


1.4 


0.06 


1.03 


45.7 


2.15 



Here ks is the Boltzmann constant. Then, the value of 
the temperature of the thermal bath can be obtained by 
multiplying the values of H(0)/kB by the reduced tem- 
perature. For instance, in the a-helix case the reduced 
temperature T — 0.061419 (see caption Fig. ||) corre- 
sponds to 0.1 K of the thermal bath if v(0) — 1.003 or 
3.8 K if u(0) = 1.03. 



V. SUMMARY AND CONCLUSIONS 

We have studied the non-equilibrium diffusion dynam- 
ics of lattice solitons on a classical chain of atoms under 
thermal fluctuations, namely soliton dynamics when the 
system energy is close to its stationary value. The inter- 
action potential between the atoms is harmonic plus a 
cubic anharmonicity. The chain is coupled to a thermal 
bath at a given temperature. For that reason we have in- 
cluded dissipation and noise in the discrete equations of 
motion of the chain. Here, it is assumed that the energy 
dissipation is provided by the irreversible processes aris- 
ing from the finite velocity of the relative displacements 
between particles in the chain. Thus the dissipative term 
takes the form of a hydrodynamical damping which is 
extensively used in elasticity theory. The noise term 
which fulfills the fluctuation-dissipation theorem then be- 
comes a discrete gradient of white noise delta-correlated 
in space and time. In the continuum approach our origi- 
nal discrete set of equations leads to a form of noisy KdV- 
Burgers equation. At this point we have used a collective 
coordinate approach to study the diffusion dynamics of 
both position and velocity of the soliton. The soliton po- 
sition and the inverse soliton width have been found to 
be good collective coordinates to describe the soliton dif- 
fusion. We have derived two stochastic ordinary differen- 
tial equations with multiplicative noise which have been 



solved analytically using stochastic perturbation analy- 
sis. 

For low-energy solitons, whose velocities are close to 
the sound velocity, our molecular dynamics simulation 
has confirmed our analytical predictions. Namely, nor- 
mal diffusion of lattice solitons governs short times, while 
superdiffusive behavior is present for long times. The 
time range of the normal diffusion depends on the initial 
velocity of the soliton: it is large for velocities close to 
the sound velocity and short for high velocities. The col- 
lective coordinate approach does not take into account 
the noise-induced phonon bath, however we have shown 
that this does not play an important role except when 
the reduced temperature (temperature in units of the 
initial soliton energy H(0)) of the thermal bath is high. 
In that regime the soliton diffusion is normal. In this 
case, for a given temperature, we have estimated in our 
simulations the value of the diffusion constant due to the 
noise-induced phonon bath when the system energy is 
close to its stationary value. We have subtracted this 
value from the full value of the diffusion constant which 
is not only due to the induced phonons but also due to the 
noise. The order of magnitude of the resultant value of 
this subtraction is predicted by the collective coordinate 
approach. 

Since we do not observe in our numerical results any 
dependence on the size of the system, we may expect 
similar results for very large systems, i.e. for N >> 1500. 

We provided an example by using an approximation of 
the Lennard- Jones potential to determine the tempera- 
ture of the thermal bath in the cases of a-helical proteins 
and crystal inertial gases. 

Finally, our results above point out the robustness of 
lattice solitons. In fact, they can exist even for higher 
values of temperature and damping constant than those 
explicitly considered in the present article. On the other 
hand, for lower values of the temperature the variances 
of the soliton position and velocity turn out to be very 
small because they scale with the temperature. So it is 
very difficult to observe them. 
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APPENDIX A: DERIVATION OF THE NOISE 
TERM 

The goal of this Appendix is to find the form of 
the noise force which would satisfy the fluctuation- 
dissipation theorem. The associated set of Langevin 
equations of the classical chain of atoms under thermal 
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fluctuations are 

dP n 

dt 
dY n 

dt 



where 



jpDamping 



Pn 

M ' 
dH 



Nc 



^Damping 



(Al) 
(A2) 



dU 



Mv I — 2— 1 — 

dt dt dt 



(A3) 



and F mse (t), which satisfies the fluctuation-dissipation 
theorem, is to determined. P n is the momentum, Y n 
denotes longitudinal displacement from its equilibrium 
position, of n — th particle with mass M and velocity 
dY n /dt. H is the Hamiltonian 



H = K + U, 



(A4) 



where 



P. 



K = T,2Jl ' U = Y,V[Y n+l -Y n ] (A5) 



and V is an arbitrary potential which depends on the 
relative displaceme nts Y n+1 — Y n . The discrete Fourier 
transform of Eqs. (Al) and (A2) read 



dPk 
dt 

dYk 
dt 



where 



f k -v%P k + F k Nmse (t) 

h 

M ' 



% = 2 (1 - cos(fe)) 



We define 



p k Nmse (t) = vmi k (t), 

where £, k (t) i s delta-correlated white noise, 
f 6 (*)&'(*')) = D(k)5(t - 



(A6) 

(A7) 
(A8) 

(A9) 



and D(k) is unknown. 

The associated Fokker-Planck equation of Eqs.(A6) in 
the Stratonovich sense takes the form 



dtp = £ -%(T, 



where 



^7fcdp PfePfe + -——dp p k 
■ 1 2^ ^- fc 



JJ«(A-A(*Mn-n(*)) 



(A10) 



(All) 



In order to determine - P(fc) we have demanded the sta- 
tionary solution of Eq.(AlC) to be the Boltzmann distri- 
bution, namely 



p = J\f cxp I — 



H 



(A12) 



where H is defined in Eq. (A4) and J\[ is the normaliza- 
tion constant. Substituting Eq.(A12) into Eq.(AlC), it is 
straightforward to see that 



D(k) = 2v%Mk B T. 



(A13) 



Therefore, from Eq. (A8) together with ( A13 ), it is easy 
to show that in position space 

(F^ mse (t)F r y° ise {t')) = 

-IvM k B T6(t - t')(S n+hn , - 2<5„,„/ + <5„_i,„/). (A14) 



Finally, the relation ( A14 ) can be satisfied by the defini- 
tion 



F 



where 



No 



\t) = y/2v Mk B T(£ n+1 (t) - Ut)), (A15) 



(Ut)^(t')) = s(t-t')s nt , 



APPENDIX B: CONTINUUM LIMIT 



(A16) 



In order to reduce Eq. (|^) to a form of noisy and 
damped KdV equation we have performed two steps. 
First, we have employed the continuum approach |^|, |2^| 
in order to obtain a form of noisy and damped Bq equa- 
tion, and then we have used the reductive perturba- 
tion technique |2^, ^9) in order to obtain the noisy and 
damped KdV equation. 



1. Noisy and damped Bq equation 

Here we have used the procedure of Pnevmatikos [Q , 
who expanded Y n ±i(t) and Y n ±2(t) in a Taylor series 
around y(x,t) 7 with x — na, where the equilibrium 
atomic spacing a is regarded as an expansion parame- 
ter. Then, collecting powers of a, Eq. (|J) together with 
Eqs. d) and (0) at °(a 4 ) takes tlie form 



8 2 y = c 2 dly+pd x yd 2 x y+hd 4 x y+va 2 d 2 x d t y+aVDd x ax,t), 
where 



,3/2 



(Bl) 
(B2) 



with properties 



(d x ax,t)) = o, 

{d x Z(x,t)d xl £(x',t')) = d x ,d x 6(x-x')6(t-t'). (B3) 
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The diffusion constant takes the form 

2vk B T 



D = 



(B4) 



Other constants are 

, Ga 
c 2 = — , 
P 

, a 3 G 



2a 2 AG 



P 



P 

p = M/a. 



(B5) 



2. Noisy and damped KdV equation 



We write Eq. (Bl) in the sound velocity frame and 
make a further approximation concerning variations in 
time. In this case we may use the reductive p ertu rbation 
technique ]39|]. Therefore we rewrite Eq. (Bl) in the 
perturbation form 

d 2 t y-c 2 d 2 x y-pd x ydly~hdty = 

k (va 2 d 2 x d t y + aVI>d x £(x, t)) , (B6) 

where we have introduced a small parameter k for con- 
venience. Afterwards, we perform the following change 
of variables 



s = na(x — ct), 
so that 

d x = nad s 

where 

P 

a 



T = n 3 (3t, u = jd s y, (B7) 

and d t — n 3 [3d T — Kacd s , (B8) 



P 



1 



7 



(B9) 



"oh 12cV6/i V6h 

We have also expressed u and £ in a perturbation series 



U = Kill + K U2 + 



(BIO) 
(Bll) 



Here the parameter n indicates the magnitude of the rate 
of change, the coefficients K and k 3 in (B7) are chosen in 
order to balance the nonlinear term, and the dispersive 
term and the time derivative are of the same order in 



k. The small noise expansion (Bll) is defined such that 
the lowest order terms of noise a nd da mpin g are of the 
same order. Substituting Eqs. flBS|), (BIO) and (Bll) 
into (|B6| ), and keeping the lowest order terms, namely 
0(k 5 ), we find that 

d T u + 6ud s u + d 3 u — v\d ss u — \fD\d s ^s, t), (B12) 
where we have set u = u% and £ = £i. 



^/6va 2 c 
Vhp 



6 a 

p 3 



(B13) 



APPENDIX C: PERTURBATION ANALYSIS 

In this Appendix we develop a perturbation approach 
to the equations (see Eqs (|39| ) and (40)) 



dS(r) = A V 2 (r)dT + e^J dWxir) 



~4V7 V ?? 3 (r) 



dr)(r) 



30^i 

30 + 7T 2 



7] 3 (T)dT 



/225(231 + 8tt 2 )L>i 
\ 112(30 + 7T 2 ) 2 



15V21 + 7T 2 
2^7(30 + 7T 2 ) 



^D 1 r 1 { T )dW 2 {T] 



(CI) 



dT 



(C2) 



We interpret Eqs. ( pjj ) and (p2j) in the Ito sense 
where the Wiener process dWi{r) — £i(r)dr with 
(£i( r )£j( T ')) = 5ij5{r — t'). We seek an asymptotic so- 
lution of the form 



S(t) = s (t) + e si(r) + ■ 
v(t) = Vo{t) +£77i(t) + • 



(C3) 



Inserting Eqs. (C3) into Eqs. (CI) and flC2j ) and collect- 
ing powers of e we get 



ds (r) = 4^(t)(Zt 
d77o(r) : 



30i^i 3 

i lo dr 



30 + 7T 2 



(C4) 
(C5) 



S// 1 .(r),/ l |-W- + ^ 1 /-^ J dn- J (r) (C(i) 



90z/i , 

~30 + 7T 2 " 

15^21 + 



4^ V 

Voi^mi^dT + 



2^(30 + 7T 2 ) 
Solving Eqs. (|cj) and ( |C5| ) we obtain 

,% 2 (o) 



«o(t) = 4 
%(t) = 



A 

Vo(0) 



log (1 + At) 



with 



A 



VTTA7 
60^(0) 

30 + 7T 2 



(C7) 

(C8) 
(C9) 

(C10) 



Inserting Eq. (|C9|) in (C7) and solving, with the initial 
condition rji(0) = 0, we find 



45 Di(231 + 8tt 2 )((1 + At) 5 / 2 - 1) 
56(30 + tt 2 ) 2 A(1 + At) 3 / 2 



15V21 + TT 2 
2^/7(30 + 7T 2 ) 



y/D^tfp) (l + \T'f/ 4 dW 2 (T'). 



(Cll) 
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FIG. 6: Var(x) (Eq. (014)) and Var(t>) (Eq.(C17)) vs. time compared with results from a numerical solution of Eqs. 
and (|38|). In panels (a) and (b) T = 5 x 10~ , and in panels (c) and (d) T = 5 x 1CP 5 . A and B in all cases corresponcPto 
v(0) = 1.005 and v(0) = 1.001, respectively, v = 0.003. Solid line: analytical prediction, dashed line: numerical solution. 



Then inserting Eqs. flCll| ) and flCJ) in Eq flCJ) and 
integrating once we get 

15^(231 + 8^h (0) / s/2 
Sl(T)= -— (1 + Ar) [ 2{1 + XT) ~ 



5At-2 



7(30 + tt 2 ) 2 A 2 (1 



4v/7% 3/2 (0) 



(1 + Ar') 3/4 dW 1 {r') + 



60V21 + 7T 2 
77(30 + 7T 2 ) 

T 

dr>- 1 



Vl + At' 



DiVo /2 (0) x 

(l + Ar") 5/4 dVF 2 (r"). (C12) 



To this order we have Var(S(r)) = e 2 Var(si(T)) and 
thus we finally obtain 



Var(S(r)) 



112r, 3 (0) 448r, 3 (0) 



2I 75 D\ _ , 225I? 1 A t2+0(t3) 

(C13) 



The full expressions of Var(S) to this order of perturba- 
tion is given by Eqs. (f42|). Notice that the variance in 
the rest frame reads 

Var(x(t)) = Xvar(S(/3t)). (C14) 
or 



Concerning soliton velocity, up to first order perturbation 
it reads 



4?7 2 (r) = 4t7 2 (t) + 8e % (t^t). 



(C15) 



Then, by substituting Eqs. (jC9|) and ( |CTl| ) in Eq. ( |C15| ) 
it is straightforward to see that 

Var(Ar] 2 (T)) = 64 e 2 {^{T)) 2 Var{rji{T)) 



2 / 3600(21 + 7T 2 )£>i 
* ^ 7(30 + 7T 2 ) 2 



% 3 (0)r 



9900(21 + tt 2 )£>i A 
7(30 + 7T 2 ) 2 



% J (0)^+O(r J ) . (C16) 



The full expressions of VartArj 2 ) to this order of per- 
turbation is given by Eqs. (|42j). In the rest frame this 
variance reads 



Var(v(t)) = ( ^- ) Var(ir) 2 (f3t)). (C17) 



In Fig. ^| we show some examples of Var(x) and Var(v) 
compared with results from a numerical solution of Eqs. 
(|37| ) and ( |38| ) for which we have used the Heun method 
p5| . The variances have been obtained by averaging over 
1000 realizations. 
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and the other one due to the soliton, (H^° p (t)). 



1000 2000 3000 4000 5000 



FIG. 7: {H^ q )/Nk B T vs. t. v = 0.003. 



APPENDIX D: THERM ALIZATION PROCESS 

From the generalized equipartition theorem pp[ we 
have that 



(Dl) 



when the system is in thermal equilibrium wit h an exter- 
nal bath at temperature T. The relation (Dl) is strictly 
satisfied in the harmonic limit of the Hamiltonian (Tfl) 
namely A — 0. For finite values of A the relation ( pi| ) 
is a rather good approximation to evaluate the temper- 
ature of our system when the relative displacements are 
sufficiently small, namely 



Y n+1 (t)-Y n (t) « 



3 

2A' 



(D2) 



Notice that A = 1 in our simulations. The condition ( |P2| ) 
can be obtained by comparing the harmonic term with 
the cubic anharmonicity in ([!]) and it is always satisfied 
in the present work. 

So, in order to examine the thermalization process, one 
can define the ensemble average 

(W*)> = (g^(*)g|)- (° 3 ) 
For finite times (H eqp (t)} possesses two contributions, 



(H eqp (t)) = (Ht h qp (t)) 



So 



(D4) 



Notice that for a very large system (N >> 1500) the 
contribution of the soliton energy, (H^° p (t)), can be 
neglected. However, in our system of 1500 sites this 
contribution is appreciable. In fact, the time evolution 
of (H eqp (t)) for different initial soliton energies presents 
different values due to the soliton contribution. Here, 
we remark that in thermal equilib rium , i.e. t — ► oo, 
these differences vanish. In Eq. (D4) {Hl qp {t)) gives 
us information about the time evolution of the tem- 
perature of the system in terms of the temperature of 
the thermal bath. The soliton contribution (H^° p (t)) 



can be evaluated numerically using Eq. (D3) when the 
soliton propagates in the presence of the damping but 
without noise. Then one can perform the numerical 
subtraction (H eqp {t)) - (H^ p (t)) to obtain (#*(*)), 
which is shown in a normalized form in Fig. M. We 
remark that we observe the same result here in systems 
with the double number of sites, namely 3000. The 
normalized thermalization process depends only on the 
damping constant which has the same value v = 0.003 in 
all our simulations. Notice that for times t'z'S/v = 1000 
there is a fast relaxation process, but for larger times, 
t^,3/v, the energy system approaches very slowly 
its stationary value. The thermal equilibrium of the 
system with the external bath corresponds to the case 
(Htl(t))/k B T=l, 



APPENDIX E: PROFILES 

In Fig. U we compare snapshots of the system at 
t = 5000 with and without noise and in the presence 
of damping. Figs. ||a and ||b correspond to both the 
kink shape (absolute displacements) and the pulse shape 
(relative displacements), respectively. Both shapes, with 
and without noise, in Fig. ^b were reconstructed from the 
shapes in Fig. @a, respectively, by using the algorithm 



one due to the coupling to the thermal bath, (Hl qp (t)), (| 
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FIG. 8: Lattice soliton profiles: pulse shape (a), kink shape (b). Solid line: noisy shape, Dashed line: shape in the presence of 
damping only. T = 5 x 10~ 5 , v(0)=1.003, v = 0.003 and t = 5000 
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